#include <math.h>

#define FLAG_DIRECT 1
#define FLAG_INVERSE -1
#define NEG_2_PI -6.283185307179586476925286766559

void DFT(int F, int L, double* X, double* Y) {
    int I, J;
    double RAD, MR, MI;
    double X0[L], Y0[L];
    for (I = 0; I < L; ++I) {
        X0[I] = X[I];
        Y0[I] = Y[I];
    }
    for (I = 0; I < L; ++I) {
        X[I] = 0;
        Y[I] = 0;
        for (J = 0; J < L; ++J) {
            RAD = NEG_2_PI * F * I * J / L;
            MR = cos(RAD);
            MI = sin(RAD);
            X[I] += MR * X0[J] - MI * Y0[J];
            Y[I] += MR * Y0[J] + MI * X0[J];
        }
    }
    if (F == FLAG_INVERSE) {
        for (int i = 0; i < L; ++i) {
            X[i] /= L;
            Y[i] /= L;
        }
    }
}


void R2FFT(int F, int L, double* X, double* Y) {
    double TEMP;
    int H, I, J, K, M, N;
    int REV[L];
    for (I = 1, H = -1; I < L; I <<= 1) ++H;
    for (I = 0; I < L; ++I) {
        J = (REV[I >> 1] >> 1) | ((I & 1) << H);
        REV[I] = J;
        if (I < J) {
            TEMP = X[I];
            X[I] = X[J];
            X[J] = TEMP;
            TEMP = Y[I];
            Y[I] = Y[J];
            Y[J] = TEMP;
        }
    }
    double RAD, WNR, WNI, WR, WI, UR, UI, TR, TI;
    for (I = 2; I <= L; I <<= 1) {
        H = I >> 1;
        RAD = NEG_2_PI * F / I;
        WNR = cos(RAD);
        WNI = sin(RAD);
        for (J = 0; J < L; J += I) {
            WR = 1;
            WI = 0;
            for (K = 0; K < H; ++K) {
                M = J + K;
                N = M + H;
                UR = X[M];
                UI = Y[M];
                TR = WR * X[N] - WI * Y[N];
                TI = WR * Y[N] + WI * X[N];
                X[M] = UR + TR;
                Y[M] = UI + TI;
                X[N] = UR - TR;
                Y[N] = UI - TI;
                TEMP = WR;
                WR = TEMP * WNR - WI * WNI;
                WI = TEMP * WNI + WI * WNR;
            }
        }
    }
    if (F == FLAG_INVERSE) {
        for (I = 0; I < L; ++I) {
            X[I] /= L;
            Y[I] /= L;
        }
    }
}

void R4FFT(int F, int L, double* X, double* Y) {
    int J, K, M, N, Q;
    int I0, I1, I2, I3;
    double RAD, RAD1, RAD2, RAD3;
    double W1R, W2R, W3R, W1I, W2I, W3I;
    double F1R, F1I, F2R, F2I, F3R, F3I, F4R, F4I;
    for (M = 1, Q = L; M < L; M <<= 2) {
        N = Q;
        RAD = NEG_2_PI * F / N;
        for (K = 0, Q >>= 2; K < Q; ++K) {
            RAD1 = K * RAD;
            RAD2 = RAD1 + RAD1;
            RAD3 = RAD1 + RAD2;
            W1R = cos(RAD1);
            W1I = sin(RAD1);
            W2R = cos(RAD2);
            W2I = sin(RAD2);
            W3R = cos(RAD3);
            W3I =sin(RAD3);
            for (I0 = K; I0 < L; I0 += N) {
                I1 = I0 + Q;
                I2 = I1 + Q;
                I3 = I2 + Q;
                F1R = X[I0] + X[I2];
                F3R = X[I0] - X[I2];
                F1I = Y[I0] + Y[I2];
                F3I = Y[I0] - Y[I2];
                F2R = X[I1] + X[I3];
                F4R = X[I1] - X[I3];
                F2I = Y[I1] + Y[I3];
                F4I = Y[I1] - Y[I3];
                X[I0] = F1R + F2R;
                F2R = F1R - F2R;
                F1R = F3R - F4I;
                F3R = F3R + F4I;
                Y[I0] = F1I + F2I;
                F2I = F1I - F2I;
                F1I = F3I + F4R;
                F3I = F3I - F4R;
                X[I2] = W2R * F2R - W2I * F2I;
                Y[I2] = W2R * F2I + W2I * F2R;
                if (F == FLAG_INVERSE) {
                    X[I1] = W1R * F1R - W1I * F1I;
                    Y[I1] = W1R * F1I + W1I * F1R;
                    X[I3] = W3R * F3R - W3I * F3I;
                    Y[I3] = W3R * F3I + W3I * F3R;
                } else {
                    X[I1] = W1R * F3R - W1I * F3I;
                    Y[I1] = W1R * F3I + W1I * F3R;
                    X[I3] = W3R * F1R - W3I * F1I;
                    Y[I3] = W3R * F1I + W3I * F1R;
                }
            }
        }
    }
    double TEMP;
    for (J = 1, K = 0, Q = L >> 2, N = L - 1; J < N; ++J) {
        M = Q;
        while (K >= 3 * M) {
            K -= 3 * M;
            M >>= 2;
        }
        K += M;
        if (J < K) {
            TEMP = X[J];
            X[J] = X[K];
            X[K] = TEMP;
            TEMP = Y[J];
            Y[J] = Y[K];
            Y[K] = TEMP;
        }
    }
    if (F == FLAG_INVERSE) {
        for (K = 0; K < L; ++K) {
            X[K] /= L;
            Y[K] /= L;
        }
    }
}
